Introduction

In this lab we’ll introduce the add-on package tmap, which can be used to produce much nicer maps than the default mapping functions in R. Before starting the lab, make sure this is installed.

Remember to use RStudio’s script editor (top left panel) to enter your commands and then copy-paste these to the console when you are ready. This will keep a copy of the commands that you can save and return to.

Before starting, remember to create a working directory (e.g. module15). Next download the files for today’s lab from the Canvas page and move them to this directory (unzip the compressed files). You’ll need the following files:

  • Climate dataset for Western North America: WNAclimate.csv
  • New York state polygon data in a shapefile: NY_Data.zip
  • A shapefile of country borders: ne_50m_admin_0_countries.zip
  • A NetCDF file of global monthly air temperature: air.mon.ltm.nc

Then start RStudio, and change your working directory from the [Session] menu > [Set working directory] > [Choose directory…]. Remember to check that R can see the files by running the list.files() command.

The tmap package

tmap works in a very similar way to ggplot2, by building a series of layers and map geometries and elements. In general, we start by using tm_shape() to identify the spatial object to be used, and add geometries are added, including filled polygons, borders and symbols. Finally, we can add legends, scale bars, etc.

library(tmap)

Example 1: Syracuse data

For the first example, we’ll make some maps using the subset of the NY8 shapefile for Syracuse (we used this in a previous module). Start by loading the packages we need (sf to read the data, and tmap to plot it):

library(tmap)
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1

Next load the NY8 data, and subset Syracuse. We’ll also make a quick map with the default mapping function as a reference:

NY8 <- st_read("NY_data/NY8_utm18.shp")
## Reading layer `NY8_utm18' from data source `/Users/u0784726/Dropbox/Data/devtools/geog5680/15 Optional/15c tmap/NY_data/NY8_utm18.shp' using driver `ESRI Shapefile'
## Simple feature collection with 281 features and 17 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 358241.9 ymin: 4649755 xmax: 480393.1 ymax: 4808545
## projected CRS:  WGS 84 / UTM zone 18N
Syracuse <- NY8[NY8$AREANAME == "Syracuse city", ]
plot(Syracuse[, "POP8"])

Now let’s do the same things with tmap. First, let’s make a simple map showing the polygon outlines using tm_borders():

tm_shape(Syracuse) + tm_borders()

The function tm_fill() will then fill these using one of the variables in the Syracuse data set (POP8). Note that this automatically adds a legend within the frame of the figure:

tm_shape(Syracuse) + tm_borders() + tm_fill("POP8")

The color scale can be changed by setting the palette argument in tm_fill(). This includes the ColorBrewer scales described in the previous module, and the algorithm used to define the intervals. For example, to use the ‘Greens’ palette with percentile breaks:

tm_shape(Syracuse) + tm_borders() + 
  tm_fill("POP8", palette = "Greens", style = "quantile")

If you want to see the full set of palettes that you can use, install the tmaptools package and run the following code (you may also need to install shinyjs):

tmaptools::palette_explorer()

Other map elements can be added. Here we add a longitude/latitude graticule with tm_graticules(), a north arrow and a line of text with the date the map was made.

tm_shape(Syracuse) + tm_graticules(col = "lightgray") + tm_borders() + 
  tm_fill("POP8", palette = "Greens", style = "quantile") + 
  tm_compass(position = c("left", "bottom")) + 
  tm_credits("2019-10-19", position = c("right", "top"))

Interactive maps

tmap also has a function to produce a simple, interactive map. This is done by changing the mode of operation from plot to view. The following code set the mode to interactive, makes a simple plot of the Syracuse data, and resets the mode to static plots. Note that these plots will be interactive in RStudio, and can also be embedded in the html files produced in R Markdown. At the end we reset the view back to static plots. Note that you can use this interactive function with any data with latitude/longitude coordinates.

tmap_mode("view")
names(Syracuse)
tm_shape(Syracuse) + tm_borders() + tm_fill("Cases", palette = "Greens")
tmap_mode("plot")

Example 2

We’ll next make some maps with the Western North American site data. We’ll first need to read the data in:

wna_climate <- read.csv("WNAclimate.csv")
str(wna_climate)
## 'data.frame':    2012 obs. of  8 variables:
##  $ WNASEQ  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ LONDD   : num  -108 -105 -103 -110 -114 ...
##  $ LATDD   : num  50.9 55.3 41.7 44.3 59.2 ...
##  $ ELEV    : int  686 369 1163 2362 880 900 1100 1480 651 725 ...
##  $ totsnopt: num  78.8 145.4 42.7 255.1 164.9 ...
##  $ annp    : int  326 499 450 489 412 451 492 606 443 455 ...
##  $ Jan_Tmp : num  -13.9 -21.3 -4.2 -10.9 -23.9 -17.5 -18.8 -11 -22.4 -23.5 ...
##  $ Jul_Tmp : num  18.8 16.8 23.3 14.1 14.4 13.8 13.2 12.8 15.3 15.4 ...

Now let’s convert this to a sf object (a points data frame):

wna_climate <- st_as_sf(wna_climate, 
                        coords = c("LONDD", "LATDD"),
                        crs = 4326)

Individual symbols can be plotted on a color scale using tm_symbols.

tm_shape(wna_climate) + tm_symbols(col="Jan_Tmp")

This takes the same arguments as tm_fill() for the color palette. We’ll use a red to blue color scale from RColorBrewer. The minus sign before the palette name reverses the order of the colors. As there is a large amount of overlap between the sites, we also add an alpha level to make the symbols transparent.

tm_shape(wna_climate) + 
  tm_symbols(col="Jan_Tmp", alpha = 0.5, palette = "-RdBu")

Next, load the Natural Earth shapefile of country boundaries.

countries <- st_read("ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp")

We can now add this to our existing plot. Note that as this is a different spatial object, we have to use tm_shape() a second time to reference this, then use tm_borders() to add the lines.

tm_shape(wna_climate) + 
  tm_symbols(col="Jan_Tmp", alpha = 0.75, palette = "-RdBu") + 
  tm_shape(countries) + tm_borders(col="gray")

We can use tm_style() to alter the appearance of the map:

tm_shape(wna_climate) + 
  tm_symbols(col="Jan_Tmp", alpha = 0.75, palette = "-RdBu", size = 0.5) + 
  tm_shape(countries) + tm_borders(col="gray") + tm_style("cobalt")

Like ggplot2, the figures that are made can be saved as R objects. Here we make two maps (January and July temperature) and save them as tm1 and tm2

tm1 <- tm_shape(wna_climate) + 
  tm_symbols(col="Jan_Tmp", alpha = 0.75, palette = "-RdBu", size = 0.5) + 
  tm_shape(countries) + tm_borders(col="gray") + tm_style("classic") 

tm2 <- tm_shape(wna_climate) + 
  tm_symbols(col="Jul_Tmp", alpha = 0.75, palette = "-RdBu", size = 0.5) + 
  tm_shape(countries) + tm_borders(col="gray") + tm_style("classic")

We can now use tm_arrange() to make a single figure with the two maps:

tmap_arrange(tm1, tm2)
## Variable(s) "Jan_Tmp" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
## Variable(s) "Jan_Tmp" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Example 3

In the final example, we’ll make figures using the global air temperature dataset. Start by reading the data using the raster package:

library(raster)
r <- raster("air.mon.ltm.nc", varname="air")
r <- rotate(r)
names(r)
## [1] "Monthly.Long.Term.Mean.Air.Temperature.at.sigma.level.0.995"

The variable name in our raster layer is quite long and complex, so to reduce the amount of typing in these examples, we’ll change ti to something simpler

names(r) <- "jan_tmp"

Now let’s make a plot using tm_raster(), which again takes similar options for color palettes. We’ll also add the country borders. Note that there is also a function tm_rgb() that allows you to plot multi-band images.

tm_shape(r) + 
  tm_raster("jan_tmp", style="fisher", palette="-RdBu") +
  tm_shape(countries) + tm_borders()

This throws an error about the projections. We’ll simply set the projection of the raster r using the crs() function:

crs(r) <- 4326

And now remake the figure

tm_shape(r) + 
  tm_raster("jan_tmp", palette="RdBu") +
  tm_shape(countries) + tm_borders()

There’s a few things wrong here. - The color palette is the wrong way around, so we can reverse this by prepending a - before the palette name. - The palette intervals are quite large, so we don’t see much of the spatial variation, particularly over mountain chains. We add n=11 to the function to increase the number of intervals - The legend doesn’t show up very well

We can address this last point in a couple of ways. First by adding a background to the legend:

tm_shape(r) + 
  tm_raster("jan_tmp", style="fisher", palette="-RdBu", n = 9, title = "January temperature") +
  tm_shape(countries) + tm_borders() + 
  tm_layout(legend.bg.color = "white", legend.bg.alpha = 0.6)

Alternatively, we can do this a little by moving the color legend outside of the plotting area. We’ll increase the number of color classes to 9, and add a histogram showing the frequency of different values.

tm_shape(r) + 
  tm_raster("jan_tmp", style="fisher", palette="-RdBu", legend.hist = TRUE, 
            n = 9, title = "January temperature") +
  tm_shape(countries) + tm_borders() + tm_layout(legend.outside = TRUE, 
                                                 legend.outside.position = "left")

## [1] "#134B86" "#3581B9" "#82BBD8" "#CAE1EE" "#F7F7F7" "#FDDBC7" "#F4A582"
## [8] "#D6604D" "#B2182B"

Exercise

There is no additional exercise for today, but you need to submit your R script or (preferably) a Markdown document with the code examples from the lab